% Scripts to read data and compute magnetic transfer functions

% script to process IPM and PPT data

% The data are 1 month 1 sec for Feb 2010
ipm = load('/Users/manojnair/projects/tsunami/IPM_PPT_2010_0910_mean_removed/IPM.nomean');
ppt = load('/Users/manojnair/projects/tsunami/IPM_PPT_2010_0910_mean_removed/PPT.nomean');
%define time axis
t = linspace(datenum(2010,09,1,0,0,0),datenum(2010,10,31,23,59,59),61*24*3600);

L = ppt > 99999.0;
ppt(L) = NaN;


L = ipm > 99999.0;
ipm(L) = NaN;

L = isnan(ipm(:,1)) | isnan(ipm(:,2)) | isnan(ipm(:,3));


ipm(:,1) = interp1(t(~L), ipm(~L,1), t, 'linear');
ipm(:,2) = interp1(t(~L), ipm(~L,2), t, 'linear');
ipm(:,3) = interp1(t(~L), ipm(~L,3), t, 'linear');


L = isnan(ppt(:,1)) | isnan(ppt(:,2)) | isnan(ppt(:,3));


ppt(:,1) = interp1(t(~L), ppt(~L,1), t, 'linear');
ppt(:,2) = interp1(t(~L), ppt(~L,2), t, 'linear');
ppt(:,3) = interp1(t(~L), ppt(~L,3), t, 'linear');


%% remove solar tidal signals


%define solar and lunar periods
periods = [4.0000    4.8000    6.0000    8.0000   11.9672   12.0000   12.4210   12.6583   23.9345   24.0000 ];
data_array = ipm(:,1);
%create a model
m=[cos(2*pi*t(:)/(periods(1)/24) ) sin(2*pi*t(:)/(periods(1)/24) )];

for j=2:length(periods)
    m = [m cos(2*pi*t(:)/(periods(j)/24) ) sin(2*pi*t(:)/(periods(j)/24) )];
end

% %fit with the data
% [s,stats] = robustfit(m, data_array);
% 
% % amplitude and phase
% P = angle(complex(s(3:2:end),s(2:2:end))); % 
% A = abs(complex(s(3:2:end),s(2:2:end)));  % BUG found. I had used arctan (x/y) 
% 
% 
% % remove unreliable tidal lines
% p_re = stats.p(2:2:end);
% p_im = stats.p(3:2:end);
% 
% L = p_re > 0.05 | p_im > 0.05;
% A(L) = 0;
% 

% LS fitting
a = m\data_array;

P = angle(complex(a(2:2:end),a(1:2:end))); % 
A = abs(complex(a(2:2:end),a(1:2:end)));  % BUG found. I had used arctan (x/y) 

% create the synthetic time series

if exist('y'),
clear y
end;

for i = 1:length(P),
    y(i,:) = A(i) * sin(2*pi*t*24/periods(i) + P(i)) ;
end;


data_corrected = data_array - sum(y); % do this for all the components
% and replace the data vector  with corrected data.

%% Bivariate TF


nfft = 4096 * 5 ;

[HxHx,   F] = pwelch(ppt(:,1),hanning(nfft),1,nfft,1);
[HyHy,   F] = pwelch(ppt(:,2),hanning(nfft),1,nfft,1);
[HxHy,   F] = cpsd(ppt(:,1),ppt(:,2),hanning(nfft),1,nfft,1);
[HzHx,   F] = cpsd(ppt(:,3),ppt(:,1),hanning(nfft),1,nfft,1);
[HzHy,   F] = cpsd(ppt(:,3),ppt(:,2),hanning(nfft),1,nfft,1);
[HyHx,   F] = cpsd(ppt(:,2),ppt(:,1),hanning(nfft),1,nfft,1);

% Txy = (HxHx .* (HzHy) - HxHy .* (HzHx))./(HxHx .* (HyHy) - HxHy .* (HyHx));
% 
% subplot(211);
% loglog((1./F), abs(Txy),'b-', 'LineWidth', 2);
% axis([ 50       10000       -inf 2]);
% 
% subplot(212);
% semilogx((1./F), unwrap(angle(Txy)) * 180/pi,'b-', 'LineWidth', 2);
% axis([ 50       10000       -180 180]);
% 
% 

for i = 1:length(HzHx);
    
EH=[HzHx(i) HzHy(i)];

HH=[HxHx(i) HxHy(i); HyHx(i) HyHy(i)];

Z(i,:) = EH/HH;

end;


%% Plot

L = isinf((1./F)) | (1./F) < 50;

subplot(211);
loglog((1./F(~L)), abs(Z(~L,2)),'b*', 'LineWidth', 1);
axis([ 50       10000       -inf 2]);
grid on
hold on
loglog((1./F(~L)), abs(Z(~L,1)),'r*', 'LineWidth', 1);
subplot(212);
semilogx((1./F(~L)), angle(Z(~L,2))* 180/pi,'b*', 'LineWidth', 1);
hold on
semilogx((1./F(~L)), angle(Z(~L,1))* 180/pi,'r*', 'LineWidth', 1);
axis([ 50       10000       -180 180]);

grid on

%% Now Plot the data

subplot(211);
set(gca, 'FontSize', 16);
%legend('TF1 (Matlab)','TF2 (Matlab)','TF1 (EMTF)','TF2 (EMTF)');
legend('TF1','TF2');

xlabel('Period (S)');
ylabel(' TF Magnitude');
subplot(212);
set(gca, 'FontSize', 16);
%legend('TF1 (Matlab)','TF2 (Matlab)','TF1 (EMTF)','TF2 (EMTF)');
legend('TF1','TF2');

xlabel('Period (S)');
ylabel(' TF Phase');
subplot(211);
%title('Bx_{ppt} = TF1*Bx_{ipm} + TF2*By_{ipm}');
%title('Bz_{IPM} = TF1*Bx_{IPM} + TF2*By_{IPM}');
%title('Bx_{IPM} = TF1*Bx_{PPT} + TF2*By_{PPT}')
%title('Bz_{PPT} = TF1*Bx_{IPM} + TF2*By_{IPM}')
%title('Bz_{IPM} = TF1*Bx_{PPT} + TF2*By_{PPT}')
title('Bz_{PPT} = TF1*Bx_{PPT} + TF2*By_{PPT}')


subplot(212);
%title('Bx_{ppt} = TF1*Bx_{ipm} + TF2*By_{ipm}');
%title('Bz_{IPM} = TF1*Bx_{IPM} + TF2*By_{IPM}');
%title('Bz_{IPM} = TF1*Bx_{PPT} + TF2*By_{PPT}')
%title('Bz_{PPT} = TF1*Bx_{IPM} + TF2*By_{IPM}')
%title('Bz_{IPM} = TF1*Bx_{PPT} + TF2*By_{PPT}')
title('Bz_{PPT} = TF1*Bx_{PPT} + TF2*By_{PPT}')



%% Now compare it with data from EMTF

cfile= '/Users/manojnair/projects/tsunami/x2y2z2z2z1.zss'

TF = read_egbert_tf(cfile);
% 1=ipm and 2= ppt

% plot the data

Z_emtf = TF.TF;
periods = TF.T;
err = TF.SIG_s;

subplot(211);

loglog(periods, abs(Z_emtf(3,:)),'r.-');

hold on;
 for i = 1: length(Z_emtf);
    plot([periods(i) periods(i)], [abs(Z_emtf(1,i)) - ... 
        abs(err(1,i))/abs(Z_emtf(1,i)) abs(Z_emtf(1,i)) + abs(err(1,i))/abs(Z_emtf(1,i))],'k','LineWidth',1);
    hold on;
end;

hold on;
loglog(periods, abs(Z_emtf(4,1:25)),'b.-');
% loglog(periods, abs(Z(3,1:25)),'k');
% loglog(periods, abs(Z(4,1:25)),'c');
% legend('TF1','TF2','TF3','TF4');

subplot(212);

semilogx(periods, angle(Z_emtf(3,1:25)) * 180/pi,'r.-');
hold on;
semilogx(periods, angle(Z_emtf(4,1:25)) * 180/pi,'b.-');
% semilogx(periods, angle(Z(3,1:25)) * 180/pi,'k');
% semilogx(periods, angle(Z(4,1:25)) * 180/pi,'c');
% legend('TF1','TF2','TF3','TF4');

